VIS and water analysis
This knitr document contains the code used to analyze and visualize data from the VIS camera system and water data. Knit-ing this document will generate an HTML document that contains both the embedded R code and the output.
Attach the required R packages
library(qvalue, warn.conflicts = FALSE)
library(ggplot2, warn.conflicts = FALSE)
library(lubridate, warn.conflicts = FALSE)
library(MASS, warn.conflicts = FALSE)
library(car, warn.conflicts = FALSE)## Warning: package 'car' was built under R version 3.1.2library(nlme, warn.conflicts = FALSE)
library(mvtnorm, warn.conflicts = FALSE)
library(grid, warn.conflicts = FALSE)Soil volume water content
Use linear regression to create a simple model for using water volume to predict soil volume water content.
Download data for soil wet/dry weight and volume water content measurements.
if (!file.exists('soil_weigth_vwc.txt')) {
download.file('http://files.figshare.com/1939954/soil_weigth_vwc.txt',
'soil_weigth_vwc.txt')
}Read the soil volume water content data
vwc.data = read.table(file="soil_weigth_vwc.txt", sep="\t", header=TRUE)Calculate the average soil dry weight
mean(vwc.data$weight_dry)## [1] 72.92138Create a linear model for water volume and volume water content. Water volumes >= 260 appear to have saturated the soil water carrying capacity, so remove them from the model.
vwc.lm = lm(vwc_wet ~ water_vol, data=vwc.data[vwc.data$water_vol < 260,])
summary(vwc.lm)##
## Call:
## lm(formula = vwc_wet ~ water_vol, data = vwc.data[vwc.data$water_vol <
## 260, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.5426 -1.6196 0.3976 1.4214 4.0771
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.032863 1.389638 -2.182 0.0391 *
## water_vol 0.233197 0.008145 28.630 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.082 on 24 degrees of freedom
## Multiple R-squared: 0.9716, Adjusted R-squared: 0.9704
## F-statistic: 819.7 on 1 and 24 DF, p-value: < 2.2e-16Plot the model
vwc.plot = ggplot(vwc.data[vwc.data$water_vol < 260,], aes(x=water_vol, y=vwc_wet)) +
geom_point(size=2) +
geom_smooth(method='lm', formula=y~x) +
scale_x_continuous("Water volume (ml)") +
scale_y_continuous("Volume water content (%)") +
theme_bw() +
theme(axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold"))
print(vwc.plot)
Predict volume water contents for the water treatment groups
treatment.water = data.frame(water_vol=c(217, 144.5, 72))
treatment.water$vwc = predict(object = vwc.lm, newdata=treatment.water)
print(treatment.water)## water_vol vwc
## 1 217.0 47.57084
## 2 144.5 30.66407
## 3 72.0 13.75731Convert VIS camera zoom units
LemnaTec VIS camera zoom units range from 1 to 6000, which correspond to 1 to 6X zoom.
zoom.lm = lm(zoom.camera ~ zoom, data=data.frame(zoom=c(1,6000),
zoom.camera=c(1,6)))
summary(zoom.lm)##
## Call:
## lm(formula = zoom.camera ~ zoom, data = data.frame(zoom = c(1,
## 6000), zoom.camera = c(1, 6)))
##
## Residuals:
## ALL 2 residuals are 0: no residual degrees of freedom!
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.9991665 NA NA NA
## zoom 0.0008335 NA NA NA
##
## Residual standard error: NaN on 0 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: NaN
## F-statistic: NaN on 1 and 0 DF, p-value: NAVIS zoom correction
In this section we define models that are used to convert area and length between camera zoom levels to a common scale.
Download data for a reference object imaged at different zoom levels.
if (!file.exists('zoom_calibration_data.txt')) {
download.file('http://files.figshare.com/1845355/zoom_calibration_data.txt',
'zoom_calibration_data.txt')
}Read zoom calibrartion data
z.data = read.table(file="zoom_calibration_data.txt", sep="\t", header=TRUE)Calculate px per cm
z.data$px_cm = z.data$length_px / z.data$length_cmConvert LemnaTec zoom units to camera zoom units
z.data$zoom.camera = predict(object = zoom.lm, newdata=z.data)Zoom correction for area
Fit a variety of regression models to relative object area by zoom level. Test whether an exponential or 2nd order polynomial model fits best (lowest AIC)
Non-linear regression (exponential)
area.coef = coef(nls(log(rel_area) ~ log(a * exp(b * zoom.camera)),
z.data, start = c(a = 1, b = 0.01)))
area.coef = data.frame(a=area.coef[1], b=area.coef[2])
area.nls = nls(rel_area ~ a * exp(b * zoom.camera),
data = z.data, start=c(a=area.coef$a, b=area.coef$b))
summary(area.nls)##
## Formula: rel_area ~ a * exp(b * zoom.camera)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a 0.37245 0.02159 17.25 <2e-16 ***
## b 0.86668 0.01504 57.61 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3306 on 32 degrees of freedom
##
## Number of iterations to convergence: 3
## Achieved convergence tolerance: 1.605e-06Non-linear regression (polynomial)
area.pol = lm(rel_area ~ zoom.camera + I(zoom.camera^2), z.data)
summary(area.pol)##
## Call:
## lm(formula = rel_area ~ zoom.camera + I(zoom.camera^2), data = z.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.8494 -0.3820 0.1176 0.3173 0.7494
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.42793 0.48203 9.186 2.34e-10 ***
## zoom.camera -4.68583 0.40917 -11.452 1.14e-12 ***
## I(zoom.camera^2) 1.64781 0.07786 21.163 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.431 on 31 degrees of freedom
## Multiple R-squared: 0.9907, Adjusted R-squared: 0.99
## F-statistic: 1643 on 2 and 31 DF, p-value: < 2.2e-16AIC
AIC(area.nls, area.pol)## df AIC
## area.nls 3 25.16270
## area.pol 4 44.12291The exponential model minimizes AIC. Plot exponential model.
nls.plot = ggplot(z.data, aes(x=zoom.camera, y=rel_area)) +
geom_point(size=2.5) +
scale_x_continuous(lim=c(0.5,4.5), "Camera zoom setting") +
scale_y_continuous(lim=c(0,16),
"Reference object relative pixel area") +
stat_smooth(data=z.data, aes(x=zoom.camera, y=rel_area),
method="nls", se=FALSE, formula=y ~ a * exp(b * x),
start=c(a=area.coef$a, b=area.coef$b)) +
theme_bw() +
theme(axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold"))nls.plot = nls.plot + labs(title='Figure 12A')
print(nls.plot)
Zoom correction for length
Fit a variety of regression models to px/cm by zoom level. Test whether an exponential or 2nd order polynomial model fits best (lowest AIC).
Non-linear regression (exponential)
len.coef = coef(nls(log(px_cm) ~ log(a * exp(b * zoom.camera)),
z.data[z.data$camera == 'VIS SV',], start = c(a = 1, b = 0.01)))
len.coef = data.frame(a=len.coef[1], b=len.coef[2])
len.nls = nls(px_cm ~ a * exp(b * zoom.camera),
data = z.data[z.data$camera == 'VIS SV',],
start=c(a=len.coef$a, b=len.coef$b))
summary(len.nls)##
## Formula: px_cm ~ a * exp(b * zoom.camera)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a 8.370121 0.150430 55.64 <2e-16 ***
## b 0.445769 0.005504 80.98 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5106 on 14 degrees of freedom
##
## Number of iterations to convergence: 3
## Achieved convergence tolerance: 2.15e-06Length zoom correction using a 2 order polynomial (px/cm given a zoom setting). Length correction only works for side-view images right now because scale in top-view images are affected by the plant lifter position in addition to zoom.
len.poly = lm(px_cm ~ zoom.camera + I(zoom.camera^2),
data=z.data[z.data$camera == 'VIS SV',])
summary(len.poly)##
## Call:
## lm(formula = px_cm ~ zoom.camera + I(zoom.camera^2), data = z.data[z.data$camera ==
## "VIS SV", ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.29376 -0.11153 0.02836 0.14223 0.21043
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.89815 0.30540 48.78 4.13e-16 ***
## zoom.camera -3.74980 0.27285 -13.74 4.04e-09 ***
## I(zoom.camera^2) 3.12322 0.05472 57.08 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1741 on 13 degrees of freedom
## Multiple R-squared: 0.9998, Adjusted R-squared: 0.9998
## F-statistic: 3.402e+04 on 2 and 13 DF, p-value: < 2.2e-16AIC
AIC(len.nls, len.poly)## df AIC
## len.nls 3 27.762784
## len.poly 4 -5.857059The polynomial model minimizes AIC. Plot polynomial model.
poly.plot = ggplot(z.data[z.data$camera == 'VIS SV',], aes(x=zoom.camera, y=px_cm)) +
geom_point(size=2.5) +
scale_x_continuous(lim=c(0.5,4.5), "Camera zoom setting") +
scale_y_continuous(lim=c(0,50),
"Reference object length scale (px/cm)") +
stat_smooth(data=z.data[z.data$camera == 'VIS SV',],
aes(x=zoom.camera, y=px_cm), method="lm",
formula=y ~ x + I(x^2)) +
theme_bw() +
theme(axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold"))poly.plot = poly.plot + labs(title='Figure 12B')
print(poly.plot)
Height modeling
In this section we measure how well PlantCV estimates plant height.
Download data manually measured plant height data set (n = 173).
if (!file.exists('height_model_data.txt')) {
download.file('http://files.figshare.com/1845361/height_model_data.txt',
'height_model_data.txt')
}Read height data.
ht.data = read.table(file="height_model_data.txt",sep="\t",header=TRUE)Convert LemnaTec zoom units to camera zoom units
ht.data$zoom.camera = predict(object = zoom.lm, newdata=ht.data)Convert pixel traits to cm with zoom correction. px/cm calibration
px_cm = predict(object = len.poly, newdata=ht.data)
ht.data$height_above_bound_cm = ht.data$height_above_bound / px_cmHeight linear model
height.model = lm(manual_height~height_above_bound_cm, ht.data)
summary(height.model)##
## Call:
## lm(formula = manual_height ~ height_above_bound_cm, data = ht.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.9229 -0.3034 -0.1263 0.3916 5.2015
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.303813 0.149776 2.028 0.0441 *
## height_above_bound_cm 0.852158 0.002789 305.528 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.103 on 171 degrees of freedom
## Multiple R-squared: 0.9982, Adjusted R-squared: 0.9982
## F-statistic: 9.335e+04 on 1 and 171 DF, p-value: < 2.2e-16Plot the height linear model.
hlm.plot = ggplot(ht.data, aes(x=height_above_bound_cm, y=manual_height)) +
geom_point(aes(color=factor(round(zoom.camera,1))),size=2.5) +
geom_smooth(method="lm", formula=y~x, color='black') +
scale_x_continuous(lim=c(0,110), breaks=c(0,25,50,75,100), "Estimated height (cm)") +
scale_y_continuous(lim=c(0,110), breaks=c(0,25,50,75,100), "Manually measured height (cm)") +
theme_bw() +
theme(legend.position=c(0.2,0.75),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Zoom")hlm.plot = hlm.plot + labs(title='Figure 3A')
print(hlm.plot)
Biomass modeling
In this section we model fresh- and dry-weight above ground biomass using image measurements.
Download data manually measured plant biomass data set (n = 41).
if (!file.exists('biomass_model_data.txt')) {
download.file('http://files.figshare.com/1845360/biomass_model_data.txt',
'biomass_model_data.txt')
}Read biomass data.
st.data = read.table(file='biomass_model_data.txt', sep="\t", header=TRUE,
stringsAsFactors=FALSE)Create out-of-frame indicator variable.
st.data$outind = NA
st.data[st.data$out_of_frame == T,]$outind = 1
st.data[st.data$out_of_frame == F,]$outind = 0Genotype indicator variable.
st.data$group = NA
st.data[st.data$genotype == 'A10',]$group = 0
st.data[st.data$genotype == 'B100',]$group = 1Fresh-weight biomass
Full model. Includes side-view area, top-view area and height.
fw.full = lm(fresh_weight ~ sv_area * tv_area * height_above_bound, st.data)Step-wise model selection with AIC.
fw.step = stepAIC(fw.full, direction="both")## Start: AIC=-50.4
## fresh_weight ~ sv_area * tv_area * height_above_bound
##
## Df Sum of Sq RSS AIC
## - sv_area:tv_area:height_above_bound 1 0.055428 8.1731 -52.122
## <none> 8.1177 -50.401
##
## Step: AIC=-52.12
## fresh_weight ~ sv_area + tv_area + height_above_bound + sv_area:tv_area +
## sv_area:height_above_bound + tv_area:height_above_bound
##
## Df Sum of Sq RSS AIC
## - tv_area:height_above_bound 1 0.024335 8.1974 -54.000
## - sv_area:tv_area 1 0.061951 8.2350 -53.812
## - sv_area:height_above_bound 1 0.108786 8.2819 -53.580
## <none> 8.1731 -52.122
## + sv_area:tv_area:height_above_bound 1 0.055428 8.1177 -50.401
##
## Step: AIC=-54
## fresh_weight ~ sv_area + tv_area + height_above_bound + sv_area:tv_area +
## sv_area:height_above_bound
##
## Df Sum of Sq RSS AIC
## - sv_area:tv_area 1 0.31153 8.5089 -54.471
## <none> 8.1974 -54.000
## - sv_area:height_above_bound 1 0.68121 8.8786 -52.727
## + tv_area:height_above_bound 1 0.02433 8.1731 -52.122
##
## Step: AIC=-54.47
## fresh_weight ~ sv_area + tv_area + height_above_bound + sv_area:height_above_bound
##
## Df Sum of Sq RSS AIC
## <none> 8.5089 -54.471
## - sv_area:height_above_bound 1 0.42578 8.9347 -54.469
## + sv_area:tv_area 1 0.31153 8.1974 -54.000
## - tv_area 1 0.54733 9.0563 -53.915
## + tv_area:height_above_bound 1 0.27391 8.2350 -53.812summary(fw.step)##
## Call:
## lm(formula = fresh_weight ~ sv_area + tv_area + height_above_bound +
## sv_area:height_above_bound, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.19711 -0.10790 -0.00974 0.08153 1.55179
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.885e-01 2.513e-01 0.750 0.4581
## sv_area 4.673e-05 4.812e-06 9.710 1.36e-11 ***
## tv_area -1.045e-05 6.864e-06 -1.522 0.1368
## height_above_bound -3.280e-02 1.451e-02 -2.260 0.0299 *
## sv_area:height_above_bound 1.166e-07 8.688e-08 1.342 0.1879
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4862 on 36 degrees of freedom
## Multiple R-squared: 0.9835, Adjusted R-squared: 0.9816
## F-statistic: 535.9 on 4 and 36 DF, p-value: < 2.2e-16AIC model
fw.aic = lm(fresh_weight ~ sv_area + tv_area + height_above_bound +
sv_area*height_above_bound, st.data)
summary(fw.aic)##
## Call:
## lm(formula = fresh_weight ~ sv_area + tv_area + height_above_bound +
## sv_area * height_above_bound, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.19711 -0.10790 -0.00974 0.08153 1.55179
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.885e-01 2.513e-01 0.750 0.4581
## sv_area 4.673e-05 4.812e-06 9.710 1.36e-11 ***
## tv_area -1.045e-05 6.864e-06 -1.522 0.1368
## height_above_bound -3.280e-02 1.451e-02 -2.260 0.0299 *
## sv_area:height_above_bound 1.166e-07 8.688e-08 1.342 0.1879
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4862 on 36 degrees of freedom
## Multiple R-squared: 0.9835, Adjusted R-squared: 0.9816
## F-statistic: 535.9 on 4 and 36 DF, p-value: < 2.2e-16The AIC model contains tv_area and height which does not have a significant coefficient, test dropping.
fw.red = lm(fresh_weight ~ sv_area, st.data)
summary(fw.red)##
## Call:
## lm(formula = fresh_weight ~ sv_area, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.17949 -0.25179 0.02035 0.23804 1.28235
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.044e-01 1.060e-01 -2.872 0.00657 **
## sv_area 3.935e-05 8.757e-07 44.931 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5003 on 39 degrees of freedom
## Multiple R-squared: 0.981, Adjusted R-squared: 0.9806
## F-statistic: 2019 on 1 and 39 DF, p-value: < 2.2e-16Goodness of fit.
anova(fw.aic, fw.red)## Analysis of Variance Table
##
## Model 1: fresh_weight ~ sv_area + tv_area + height_above_bound + sv_area *
## height_above_bound
## Model 2: fresh_weight ~ sv_area
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 36 8.5089
## 2 39 9.7630 -3 -1.2541 1.7686 0.1706SV area model.
sv.model = lm(fresh_weight ~ sv_area, st.data)
summary(sv.model)##
## Call:
## lm(formula = fresh_weight ~ sv_area, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.17949 -0.25179 0.02035 0.23804 1.28235
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.044e-01 1.060e-01 -2.872 0.00657 **
## sv_area 3.935e-05 8.757e-07 44.931 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5003 on 39 degrees of freedom
## Multiple R-squared: 0.981, Adjusted R-squared: 0.9806
## F-statistic: 2019 on 1 and 39 DF, p-value: < 2.2e-16Plot SV model
sv.model.plot = ggplot(st.data,aes(x=sv_area/1e5, y=fresh_weight)) +
geom_smooth(method="lm", color="black", formula = y ~ x) +
geom_point(size=2.5) +
scale_x_continuous("Shoot and leaf area (x10^5 px)") +
scale_y_continuous("Fresh-weight biomass (g)") +
theme_bw() +
theme(axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold"))sv.model.plot = sv.model.plot + labs(title='Figure 4A')
print(sv.model.plot)
SV area model with out-of-frame.
sv.ind.model = lm(fresh_weight ~ sv_area + outind, st.data)
anova(sv.model, sv.ind.model)## Analysis of Variance Table
##
## Model 1: fresh_weight ~ sv_area
## Model 2: fresh_weight ~ sv_area + outind
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 39 9.7630
## 2 38 7.5498 1 2.2132 11.139 0.0019 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(sv.ind.model)##
## Call:
## lm(formula = fresh_weight ~ sv_area + outind, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.75671 -0.21287 -0.00372 0.22304 1.17010
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.499e-01 9.583e-02 -2.607 0.0130 *
## sv_area 4.028e-05 8.288e-07 48.599 <2e-16 ***
## outind -5.963e-01 1.787e-01 -3.338 0.0019 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4457 on 38 degrees of freedom
## Multiple R-squared: 0.9853, Adjusted R-squared: 0.9846
## F-statistic: 1277 on 2 and 38 DF, p-value: < 2.2e-16Plot SV model with out-of-frame.
sv.model.out.plot = ggplot(st.data, aes(x=sv_area/1e5, y=fresh_weight,
group=out_of_frame, color=out_of_frame)) +
geom_point(size=2.5) +
geom_smooth(method="lm", color="black") +
scale_x_continuous("Shoot and leaf area (x10^5 px)") +
scale_y_continuous("Fresh-weight biomass (g)") +
theme_bw() +
theme(legend.position=c(0.2,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Out-of-frame")
print(sv.model.out.plot)
SV area model with genotype
sv.gt.model = lm(fresh_weight ~ sv_area + group, st.data)
summary(sv.gt.model)##
## Call:
## lm(formula = fresh_weight ~ sv_area + group, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.32109 -0.29602 0.03093 0.21960 1.15107
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.510e-01 1.342e-01 -3.362 0.00178 **
## sv_area 3.956e-05 8.637e-07 45.802 < 2e-16 ***
## group 2.647e-01 1.542e-01 1.717 0.09419 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4883 on 38 degrees of freedom
## Multiple R-squared: 0.9824, Adjusted R-squared: 0.9815
## F-statistic: 1061 on 2 and 38 DF, p-value: < 2.2e-16Plot SV model with genotype.
sv.model.gen.plot = ggplot(st.data, aes(x=sv_area/1e5, y=fresh_weight,
group=genotype, color=genotype)) +
geom_point(size=2.5) +
geom_smooth(method="lm", color="black") +
scale_x_continuous("Shoot and leaf area (x10^5 px)") +
scale_y_continuous("Fresh-weight biomass (g)") +
theme_bw() +
theme(legend.position=c(0.2,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Genotype")
print(sv.model.gen.plot)
Dry-weight biomass
SV area model
dry.sv.model = lm(dry_weight ~ sv_area, st.data)
summary(dry.sv.model)##
## Call:
## lm(formula = dry_weight ~ sv_area, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.183938 -0.047916 0.009925 0.039315 0.195700
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.295e-02 1.479e-02 -3.58 0.000938 ***
## sv_area 4.634e-06 1.222e-07 37.93 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06981 on 39 degrees of freedom
## Multiple R-squared: 0.9736, Adjusted R-squared: 0.9729
## F-statistic: 1438 on 1 and 39 DF, p-value: < 2.2e-16Plot dry-weight side-view model
dry.sv.model.plot = ggplot(st.data,aes(x=sv_area/1e5, y=dry_weight)) +
geom_smooth(method="lm", color="black", formula = y ~ x) +
geom_point(size=2.5) +
scale_x_continuous("Shoot and leaf area (x10^5 px)") +
scale_y_continuous("Dry-weight biomass (g)") +
theme_bw() +
theme(axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold"))dry.sv.model.plot = dry.sv.model.plot + labs(title='Supplemental Figure S2')
print(dry.sv.model.plot)
SV area model with out-of-frame.
dry.sv.ind.model = lm(dry_weight ~ sv_area + outind, st.data)
summary(dry.sv.ind.model)##
## Call:
## lm(formula = dry_weight ~ sv_area + outind, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.133525 -0.027886 0.006673 0.024790 0.131393
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.393e-02 1.255e-02 -3.501 0.001202 **
## sv_area 4.789e-06 1.085e-07 44.124 < 2e-16 ***
## outind -9.869e-02 2.340e-02 -4.219 0.000147 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05836 on 38 degrees of freedom
## Multiple R-squared: 0.982, Adjusted R-squared: 0.9811
## F-statistic: 1038 on 2 and 38 DF, p-value: < 2.2e-16dry.sv.model.out.plot = ggplot(st.data, aes(x=sv_area/1e5, y=dry_weight,
group=out_of_frame, color=out_of_frame)) +
geom_point(size=2.5) +
geom_smooth(method="lm", color="black") +
scale_x_continuous("Shoot and leaf area (x10^5 px)") +
scale_y_continuous("Dry-weight biomass (g)") +
theme_bw() +
theme(legend.position=c(0.2,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Out-of-frame")
print(dry.sv.model.out.plot)
SV area model with genotypes
dry.sv.gt.model = lm(dry_weight ~ sv_area + group, st.data)
summary(dry.sv.gt.model)##
## Call:
## lm(formula = dry_weight ~ sv_area + group, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.183197 -0.049205 0.008654 0.040345 0.196263
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.413e-02 1.943e-02 -2.786 0.00828 **
## sv_area 4.636e-06 1.251e-07 37.061 < 2e-16 ***
## group 2.130e-03 2.233e-02 0.095 0.92451
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07071 on 38 degrees of freedom
## Multiple R-squared: 0.9736, Adjusted R-squared: 0.9722
## F-statistic: 701 on 2 and 38 DF, p-value: < 2.2e-16dry.sv.model.gen.plot = ggplot(st.data, aes(x=sv_area/1e5, y=dry_weight,
group=genotype, color=genotype)) +
geom_point(size=2.5) +
geom_smooth(method="lm", color="black") +
scale_x_continuous("Shoot and leaf area (x10^5 px)") +
scale_y_continuous("Dry-weight biomass (g)") +
theme_bw() +
theme(legend.position=c(0.2,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Genotype")
print(dry.sv.model.gen.plot)
Analyze VIS data
In this section we analyze geometric traits from the complete VIS data set.
Download the complete VIS data set. There were a total of 6,399 snapshots with VIS image data, but the download only includes the 6,207 snapshots that were successfully processed by PlantCV. Failed snapshots generally are those that lack a plant (empty pot controls, dead plants, etc.)
if (!file.exists('vis_snapshots.txt')) {
download.file('http://files.figshare.com/1845362/vis_snapshots.txt',
'vis_snapshots.txt')
}Read data and format for analysis
vis.data = read.table(file="vis_snapshots.txt", sep=",", header=TRUE)Planting date
planting_date = as.POSIXct("2013-11-26")Add water treatment column coded in barcodes.
vis.data$treatment <- NA
vis.data$treatment[grep("AA", vis.data$plant_id)] <- 100
vis.data$treatment[grep("AB", vis.data$plant_id)] <- 0
vis.data$treatment[grep("AC", vis.data$plant_id)] <- 16
vis.data$treatment[grep("AD", vis.data$plant_id)] <- 33
vis.data$treatment[grep("AE", vis.data$plant_id)] <- 66Add plant genotype column coded in barcodes.
vis.data$genotype <- NA
vis.data$genotype[grep("p1", vis.data$plant_id)] <- 'A10'
vis.data$genotype[grep("p2", vis.data$plant_id)] <- 'B100'
vis.data$genotype[grep("r1", vis.data$plant_id)] <- 'R20'
vis.data$genotype[grep("r2", vis.data$plant_id)] <- 'R70'
vis.data$genotype[grep("r3", vis.data$plant_id)] <- 'R98'
vis.data$genotype[grep("r4", vis.data$plant_id)] <- 'R102'
vis.data$genotype[grep("r5", vis.data$plant_id)] <- 'R128'
vis.data$genotype[grep("r6", vis.data$plant_id)] <- 'R133'
vis.data$genotype[grep("r7", vis.data$plant_id)] <- 'R161'
vis.data$genotype[grep("r8", vis.data$plant_id)] <- 'R187'Add genotype x treatment group column
vis.data$group = paste(vis.data$genotype,'-',vis.data$treatment,sep='')Remove plants that were sampled for biomass measurements.
vis.data = vis.data[!vis.data$plant_id %in% st.data$barcode,]Add calendar-time data column using the Unix-time data
vis.data$date = as.POSIXct(vis.data$datetime, origin = "1970-01-01")Calculate days after planting from planting data
vis.data$dap = as.numeric(vis.data$date - planting_date)Predict fresh- and dry-weight biomass from linear models
vis.data$fw_biomass = predict.lm(object = sv.model, newdata=vis.data)
vis.data$dw_biomass = predict.lm(object = dry.sv.model, newdata=vis.data)Plant height analysis
Plot height for S. viridis and S. italica water treatments 100% and 33% full-capacity.
height.plot = ggplot(vis.data[(vis.data$genotype == 'A10' |
vis.data$genotype == 'B100') &
(vis.data$treatment == 100 |
vis.data$treatment == 33),],
aes(x=dap, y=height_above_bound, color=factor(group))) +
geom_point(size=2.5) +
geom_smooth(method="loess",size=1) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(0,120),
breaks=c(0,20,40,60,80,100,120),
name="Estimated height (cm)") +
theme_bw() +
theme(legend.position=c(0.25,0.75),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Genotype-treatment")height.plot = height.plot + labs(title='Figure 3B')
print(height.plot)
Plot height for all genotypes.
height.facet.plot = ggplot(vis.data[vis.data$treatment == 100 |
vis.data$treatment == 33,],
aes(x=dap, y=height_above_bound, color=factor(treatment))) +
geom_point(size=1) +
geom_smooth(method="loess",size=1) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(0,120),
breaks=c(0,20,40,60,80,100,120),
name="Estimated height (cm)") +
theme_bw() +
theme(legend.position=c(0.8,0.1),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Treatment") +
facet_wrap(~ genotype, ncol = 3)height.facet.plot = height.facet.plot + labs(title='Supplemental Figure S1')
print(height.facet.plot)
Genotype differences in height
For each day, calculate the mean and 95% confidence intervals for each genotype in the control water group.
height_per_day = function(vis.data) {
dap = c()
genotype = c()
int.low = c()
int.up = c()
est = c()
# Loop through each day and genotype
for(d in min(as.integer(vis.data$dap)):max(as.integer(vis.data$dap))) {
if(d != 24) {
for(g in levels(factor(vis.data$genotype))) {
h.data = vis.data[vis.data$genotype == g & vis.data$treatment == 100 &
as.integer(vis.data$dap) == d,]$height_above_bound
if(sd(h.data) != 0) {
dap = c(dap,d)
genotype = c(genotype,g)
test = t.test(h.data)
int.low = c(int.low,test$conf.int[1])
int.up = c(int.up,test$conf.int[2])
est = c(est,test$estimate)
}
}
}
}
#drought_resp = drought_resp[-14,]
results = data.frame(dap=as.numeric(dap),
genotype=genotype,
conf.int.low=int.low,
conf.int.up=int.up,
mean=est)
return(results)
}Height per day 95% CI. Write results to file.
height.perday.results = height_per_day(vis.data)
write.csv(height.perday.results,file="height.perday.100water.csv")
print(height.perday.results)## dap genotype conf.int.low conf.int.up mean
## 1 11 A10 7.8384656 9.927418 8.882942
## 2 11 B100 5.3761414 6.677522 6.026831
## 3 11 R102 2.6078989 6.747923 4.677911
## 4 11 R128 3.4255116 6.098742 4.762127
## 5 11 R133 5.8106341 11.190037 8.500336
## 6 11 R161 5.9550014 6.742792 6.348897
## 7 11 R187 4.6001185 8.238236 6.419177
## 8 11 R20 8.7685119 12.487009 10.627761
## 9 11 R70 6.9589765 8.860251 7.909614
## 10 11 R98 11.3958077 13.499331 12.447569
## 11 12 A10 10.5706419 14.040344 12.305493
## 12 12 B100 7.1779957 9.543168 8.360582
## 13 12 R102 4.1386602 10.394107 7.266383
## 14 12 R128 4.7644938 9.358302 7.061398
## 15 12 R133 0.3769178 32.979175 16.678046
## 16 12 R161 5.5667113 10.395653 7.981182
## 17 12 R187 8.9117000 14.337900 11.624800
## 18 12 R20 -24.8325684 57.616115 16.391774
## 19 12 R70 6.2768196 11.896208 9.086514
## 20 12 R98 14.4614239 19.097635 16.779529
## 21 13 A10 14.2175454 16.845714 15.531630
## 22 13 B100 10.9868413 12.771688 11.879265
## 23 13 R102 8.1006537 12.638227 10.369440
## 24 13 R128 10.0067593 13.083564 11.545162
## 25 13 R133 15.6105089 20.434429 18.022469
## 26 13 R161 11.5326626 14.465157 12.998910
## 27 13 R187 11.0472672 15.088387 13.067827
## 28 13 R20 14.8393242 21.736455 18.287890
## 29 13 R70 12.7527374 14.743597 13.748167
## 30 13 R98 17.9034635 24.295281 21.099372
## 31 14 A10 16.0353867 20.806689 18.421038
## 32 14 B100 14.8758324 16.482712 15.679272
## 33 14 R102 9.5848380 17.897359 13.741099
## 34 14 R128 11.5441914 18.807803 15.175997
## 35 14 R133 24.4394863 26.317757 25.378622
## 36 14 R161 11.3351757 16.721924 14.028550
## 37 14 R187 15.9765377 22.156071 19.066305
## 38 14 R20 9.1279616 35.869895 22.498929
## 39 14 R70 12.8683843 21.043290 16.955837
## 40 14 R98 22.7445324 26.529747 24.637140
## 41 15 A10 20.7800983 23.377751 22.078924
## 42 15 B100 17.7777752 20.518164 19.147970
## 43 15 R102 11.0573174 17.922030 14.489674
## 44 15 R128 16.7501844 21.295836 19.023010
## 45 15 R133 27.0393189 33.948830 30.494074
## 46 15 R161 17.1755797 20.612348 18.893964
## 47 15 R187 15.8387658 22.954208 19.396487
## 48 15 R20 25.1788370 33.253535 29.216186
## 49 15 R70 18.9601854 22.863978 20.912081
## 50 15 R98 29.0498999 31.955600 30.502750
## 51 16 A10 20.0577561 25.806257 22.932007
## 52 16 B100 20.0413990 23.203604 21.622501
## 53 16 R102 14.8670912 21.144546 18.005818
## 54 16 R128 16.5582461 30.333804 23.446025
## 55 16 R133 33.0362095 35.768486 34.402348
## 56 16 R161 16.4543384 23.558591 20.006465
## 57 16 R187 22.3958994 28.061383 25.228641
## 58 16 R20 9.0056082 50.912573 29.959091
## 59 16 R70 20.2903256 31.959948 26.125137
## 60 16 R98 32.7722175 35.627625 34.199921
## 61 17 A10 27.4881214 32.177027 29.832574
## 62 17 B100 25.7375286 29.145299 27.441414
## 63 17 R102 9.3832005 41.668709 25.525955
## 64 17 R128 26.8195976 35.438632 31.129115
## 65 17 R133 35.0134141 35.046324 35.029869
## 66 17 R161 18.5118930 41.274711 29.893302
## 67 17 R187 26.7251684 31.553360 29.139264
## 68 17 R20 34.5958844 35.180457 34.888171
## 69 17 R70 27.8333370 32.092941 29.963139
## 70 17 R98 34.9761994 35.063296 35.019748
## 71 18 A10 25.4649891 34.113739 29.789364
## 72 18 B100 29.5995162 33.060621 31.330069
## 73 18 R102 23.9458634 35.780013 29.862938
## 74 18 R128 28.8943167 33.278893 31.086605
## 75 18 R161 22.7950766 33.274113 28.034595
## 76 18 R187 28.6642353 38.039275 33.351755
## 77 18 R20 34.8319033 35.217714 35.024809
## 78 18 R70 32.5635770 36.332210 34.447894
## 79 18 R98 35.0377193 35.072626 35.055172
## 80 19 A10 29.9565044 33.715470 31.835987
## 81 19 B100 28.7587854 32.088308 30.423547
## 82 19 R102 18.6004422 36.868409 27.734426
## 83 19 R128 35.0349633 35.062368 35.048666
## 84 19 R133 35.0369049 35.055222 35.046063
## 85 19 R161 27.5528428 35.523189 31.538016
## 86 19 R187 30.2853896 36.052235 33.168812
## 87 19 R20 35.0431394 35.064445 35.053792
## 88 19 R70 32.8273407 35.931809 34.379575
## 89 20 A10 30.7762330 35.328096 33.052164
## 90 20 B100 31.4095265 35.224147 33.316837
## 91 20 R102 29.8675676 36.244858 33.056213
## 92 20 R128 31.8182665 36.308301 34.063284
## 93 20 R133 35.0292025 35.062924 35.046063
## 94 20 R161 27.6114760 34.513152 31.062314
## 95 20 R187 33.7013663 35.619516 34.660441
## 96 20 R20 35.0065633 35.093660 35.050112
## 97 20 R70 31.1659064 36.873618 34.019762
## 98 20 R98 35.0336568 35.066567 35.050112
## 99 21 A10 34.7350643 43.619676 39.177370
## 100 21 B100 33.9185376 40.511408 37.214973
## 101 21 R102 20.4335385 53.700974 37.067256
## 102 21 R128 45.7878729 59.113583 52.450728
## 103 21 R133 69.5748771 75.552885 72.563881
## 104 21 R161 30.8383653 46.826065 38.832215
## 105 21 R187 41.2621776 55.259613 48.260895
## 106 21 R20 56.3428640 68.941224 62.642044
## 107 21 R70 45.2148076 52.242602 48.728705
## 108 21 R98 56.5053136 69.874481 63.189898
## 109 22 A10 42.1333456 51.279316 46.706331
## 110 22 B100 35.1894676 46.590960 40.890214
## 111 22 R102 32.4270279 54.660958 43.543993
## 112 22 R128 42.7607416 57.866538 50.313640
## 113 22 R133 64.7341342 95.014087 79.874111
## 114 22 R161 31.9080496 39.897530 35.902790
## 115 22 R187 49.2083816 68.125731 58.667056
## 116 22 R20 45.5047607 113.876490 79.690625
## 117 22 R70 37.0178551 70.222627 53.620241
## 118 22 R98 64.6229237 77.916294 71.269609
## 119 23 A10 46.1205983 57.372009 51.746303
## 120 23 B100 39.0566070 49.213411 44.135009
## 121 23 R102 22.1355965 62.113193 42.124395
## 122 23 R128 56.6189532 68.402234 62.510593
## 123 23 R133 80.1898844 87.496497 83.843190
## 124 23 R161 39.5976979 48.909886 44.253792
## 125 23 R187 44.9050613 74.178923 59.541992
## 126 23 R20 65.9802051 85.598085 75.789145
## 127 23 R70 58.6470448 67.651642 63.149343
## 128 23 R98 60.6980969 79.556520 70.127309
## 129 25 A10 55.9173941 66.531485 61.224439
## 130 25 B100 42.1436759 60.516430 51.330053
## 131 25 R102 37.5833391 74.825638 56.204488
## 132 25 R128 59.0014536 77.608288 68.304871
## 133 25 R133 91.5153533 106.509863 99.012608
## 134 25 R161 41.7599259 51.527936 46.643931
## 135 25 R187 69.4178315 80.320776 74.869304
## 136 25 R20 47.3113464 137.062553 92.186950
## 137 25 R70 54.8053892 93.180449 73.992919
## 138 25 R98 79.0174138 95.332386 87.174900
## 139 26 A10 60.7309054 74.484346 67.607626
## 140 26 B100 51.9573893 67.031311 59.494350
## 141 26 R102 29.0283221 90.704621 59.866472
## 142 26 R128 69.7432021 85.949025 77.846114
## 143 26 R133 102.1131687 102.209599 102.161384
## 144 26 R161 44.7659963 77.522223 61.144110
## 145 26 R187 64.2493428 83.643787 73.946565
## 146 26 R20 95.4818290 103.551591 99.516710
## 147 26 R70 80.8950037 90.615600 85.755302
## 148 26 R98 67.3596660 96.850149 82.104907
## 149 27 A10 68.1673027 79.966027 74.066665
## 150 27 B100 52.3930899 73.531038 62.962064
## 151 27 R102 46.4307081 91.299260 68.864984
## 152 27 R128 69.7163940 92.330224 81.023309
## 153 27 R133 102.0534520 102.623594 102.338523
## 154 27 R161 47.0038080 59.572263 53.288036
## 155 27 R187 79.9731196 87.133834 83.553477
## 156 27 R20 94.6611877 107.250057 100.955622
## 157 27 R70 67.7411474 102.378830 85.059989
## 158 27 R98 92.2825030 102.888066 97.585284
## 159 28 A10 72.4019654 84.816088 78.609027
## 160 28 B100 60.2668834 76.201220 68.234052
## 161 28 R102 38.3458163 100.670481 69.508149
## 162 28 R128 88.7635776 97.020262 92.891920
## 163 28 R133 102.0427075 102.412501 102.227604
## 164 28 R161 54.3065643 86.175685 70.241125
## 165 28 R187 68.2614296 91.127546 79.694488
## 166 28 R20 101.9335479 102.697144 102.315346
## 167 28 R70 96.9145446 101.216624 99.065584
## 168 28 R98 81.5682665 106.133203 93.850735
## 169 29 A10 77.9931551 87.733707 82.863431
## 170 29 B100 64.8732783 90.167093 77.520186
## 171 29 R102 51.4886408 103.817301 77.652971
## 172 29 R128 78.2635038 103.522285 90.892894
## 173 29 R133 102.1483036 102.853222 102.500763
## 174 29 R161 60.8247415 72.598147 66.711444
## 175 29 R187 85.3621398 92.131055 88.746598
## 176 29 R20 102.1549445 102.730696 102.442820
## 177 29 R70 81.5101446 109.747356 95.628750
## 178 29 R98 100.9210632 102.477379 101.699221
## 179 30 A10 79.7975021 93.165534 86.481518
## 180 30 B100 69.9817638 86.831826 78.406795
## 181 30 R102 52.6392582 105.567685 79.103472
## 182 30 R128 100.2996586 102.191014 101.245336
## 183 30 R133 102.1362130 102.517656 102.326935
## 184 30 R161 72.2477821 89.354608 80.801195
## 185 30 R187 75.6166966 95.754845 85.685771
## 186 30 R20 102.3143461 102.675591 102.494969
## 187 30 R70 102.0282615 102.642163 102.335212
## 188 30 R98 93.7598548 105.298398 99.529126
## 189 31 A10 87.1296524 94.989746 91.059699
## 190 31 B100 76.9270063 98.740991 87.833999
## 191 31 R102 62.8067625 104.995504 83.901133
## 192 31 R128 89.1628356 105.911162 97.536999
## 193 31 R133 102.3840680 102.756520 102.570294
## 194 31 R161 68.6018081 80.383543 74.492676
## 195 31 R187 89.8228715 96.405196 93.114034
## 196 31 R20 101.8550011 103.378296 102.616648
## 197 31 R70 87.6302155 109.560625 98.595420
## 198 31 R98 102.2831228 102.641146 102.462134
## 199 32 A10 86.6320427 101.920801 94.276422
## 200 32 B100 82.8496624 97.181019 90.015341
## 201 32 R102 66.9055230 110.225409 88.565466
## 202 32 R128 105.0325661 107.488859 106.260713
## 203 32 R133 107.1520311 107.382340 107.267186
## 204 32 R161 77.4568438 96.766731 87.111787
## 205 32 R187 85.8563435 97.724316 91.790330
## 206 32 R20 107.2913801 107.379954 107.335667
## 207 32 R70 107.0598209 107.278888 107.169355
## 208 32 R98 105.1210435 108.008120 106.564582
## 209 33 A10 96.0947018 103.735921 99.915311
## 210 33 B100 89.2731829 105.354814 97.313999
## 211 33 R102 75.8742679 110.271342 93.072805
## 212 33 R128 96.1577692 110.345571 103.251670
## 213 33 R133 107.1725115 107.361860 107.267186
## 214 33 R161 81.8858866 88.986237 85.436062
## 215 33 R187 93.8716625 99.386701 96.629182
## 216 33 R20 107.0615755 107.597308 107.329442
## 217 33 R70 92.3612921 114.279011 103.320152
## 218 33 R98 107.1300106 107.383609 107.256810Treatment differences in height
Statistical analysis for height differences. Use pairwise comparisons on S. viridis and S. italica over two-day intervals.
analyze_height = function(vis.data, genotype) {
days = c()
diff.low = c()
diff.up = c()
pvals = c()
for(day in levels(factor(as.integer(vis.data$dap)))) {
day = as.integer(day)
control = vis.data[(as.integer(vis.data$dap) == day |
as.integer(vis.data$dap) == day + 1) &
vis.data$genotype == genotype &
vis.data$treatment == 100,]$height_above_bound
drought = vis.data[(as.integer(vis.data$dap) == day |
as.integer(vis.data$dap) == day + 1) &
vis.data$genotype == genotype &
vis.data$treatment == 33,]$height_above_bound
test = t.test(x=control, y=drought)
days = c(days, day)
diff.low = c(diff.low, test$conf.int[1])
diff.up = c(diff.up, test$conf.int[2])
pvals = c(pvals, test$p.value)
}
results = data.frame(dap=as.numeric(days),
conf.int.low=diff.low,
conf.int.up=diff.up,
pvalue=pvals)
return(results)
}Test for treatment effect on two-day intervals
a10.height.results = analyze_height(vis.data, 'A10')
b100.height.results = analyze_height(vis.data, 'B100')Control for multiple testing by controlling the FDR
qvalues.height = qvalue(c(a10.height.results$pvalue, b100.height.results$pvalue))
a10.height.results$qvalue = qvalues.height$qvalues[1:nrow(a10.height.results)]
b100.height.results$qvalue = qvalues.height$qvalues[
(nrow(a10.height.results)+1):(nrow(a10.height.results) + nrow(b100.height.results))]Assign genotypes for merged table
a10.height.results$genotype = 'A10'
b100.height.results$genotype = 'B100'
print(a10.height.results)## dap conf.int.low conf.int.up pvalue qvalue genotype
## 1 11 1.3757087 5.139721 1.622327e-03 1.057500e-03 A10
## 2 12 -0.7380484 4.515406 1.490756e-01 5.492423e-02 A10
## 3 13 -0.5912495 3.543282 1.562446e-01 5.516696e-02 A10
## 4 14 -1.1319577 4.329657 2.381965e-01 6.960212e-02 A10
## 5 15 -0.1299653 5.167933 6.143652e-02 3.062407e-02 A10
## 6 16 -1.3607720 7.080051 1.715524e-01 5.519138e-02 A10
## 7 17 -2.9248819 5.400151 5.353411e-01 1.191457e-01 A10
## 8 18 -0.4959663 6.673738 8.769306e-02 3.911081e-02 A10
## 9 19 1.4409925 7.969806 6.845578e-03 4.143499e-03 A10
## 10 20 3.6149867 12.029236 6.790664e-04 4.795304e-04 A10
## 11 21 6.7803624 17.104831 5.343934e-05 6.694424e-05 A10
## 12 22 10.2910019 21.102986 1.471517e-06 1.246954e-05 A10
## 13 23 9.2657053 24.761796 2.155551e-04 1.660545e-04 A10
## 14 25 7.2790852 20.484417 1.699507e-04 1.600168e-04 A10
## 15 26 8.5545706 23.844353 1.911306e-04 1.619628e-04 A10
## 16 27 11.5264141 26.176502 1.946187e-05 3.758465e-05 A10
## 17 28 10.7457158 25.976630 5.606520e-05 6.694424e-05 A10
## 18 29 11.4473342 24.918441 7.632409e-06 3.233827e-05 A10
## 19 30 11.9787099 27.114832 2.217662e-05 3.758465e-05 A10
## 20 31 12.9567838 28.606319 1.247030e-05 3.522418e-05 A10
## 21 32 12.2483513 29.889421 6.320013e-05 6.694424e-05 A10
## 22 33 3.4267550 38.584764 2.690657e-02 1.518525e-02 A10print(b100.height.results)## dap conf.int.low conf.int.up pvalue qvalue genotype
## 1 11 -2.1888257 0.4486175 0.18236609 0.05519138 B100
## 2 12 -2.2019131 1.9489490 0.89939162 0.17724157 B100
## 3 13 -3.7808622 0.5052064 0.12627973 0.05095650 B100
## 4 14 -2.7496601 2.5699624 0.94357005 0.18172165 B100
## 5 15 -4.8854099 0.8999123 0.16277335 0.05517324 B100
## 6 16 -2.8974999 3.9598806 0.74490405 0.15029213 B100
## 7 17 -5.0026064 2.9536026 0.58132070 0.12014816 B100
## 8 18 -1.9644558 5.1388196 0.35534030 0.09713326 B100
## 9 19 -2.3699180 4.9749886 0.46324842 0.10904271 B100
## 10 20 -0.8187179 6.7951985 0.11836288 0.05014997 B100
## 11 21 -8.9817009 4.6032442 0.50430459 0.11549851 B100
## 12 22 -8.3643015 4.7370840 0.57202885 0.12014816 B100
## 13 23 -11.0131272 5.4448269 0.46303807 0.10904271 B100
## 14 25 -13.4045213 2.5892801 0.17774764 0.05519138 B100
## 15 26 -17.3140056 -1.0219101 0.02867192 0.01518525 B100
## 16 27 -14.3161656 2.1808693 0.14381616 0.05492423 B100
## 17 28 -17.0011012 1.2317531 0.08754828 0.03911081 B100
## 18 29 -11.9653330 5.2377910 0.43167803 0.10904271 B100
## 19 30 -12.3111302 4.5263257 0.35381611 0.09713326 B100
## 20 31 -4.5854432 10.1681899 0.44686134 0.10904271 B100
## 21 32 -5.5221448 10.1865601 0.54834969 0.11914565 B100
## 22 33 -6.1404245 14.8562724 0.38830762 0.10282794 B100Output the A10 and B100 tables to the same file
write.table(a10.height.results, file='height.stats.csv', sep = ',',
row.names = FALSE, append = FALSE)
write.table(b100.height.results, file='height.stats.csv', sep = ',',
row.names = FALSE, append = TRUE, col.names = FALSE)Biomass analysis
Plot fresh-weight biomass for S. viridis and S. italica water treatments 100% and 33% full-capacity.
biomass.plot = ggplot(vis.data[(vis.data$genotype == 'A10' |
vis.data$genotype == 'B100') &
(vis.data$treatment == 100 |
vis.data$treatment == 33),],
aes(x=dap, y=fw_biomass, color=factor(group))) +
geom_point(size=2.5) +
geom_smooth(method="loess",size=1) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(-1,41),
name="Estimated biomass (g)") +
theme_bw() +
theme(legend.position=c(0.25,0.75),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Genotype-treatment")biomass.plot = biomass.plot + labs(title='Figure 4B')
print(biomass.plot)
Plot fresh-weight biomass for all genotypes.
biomass.facet.plot = ggplot(vis.data[vis.data$treatment == 100 |
vis.data$treatment == 33,],
aes(x=dap, y=fw_biomass, color=factor(treatment))) +
geom_point(size=1) +
geom_smooth(method="loess",size=1) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(-1,41),
name="Estimated biomass (g)") +
theme_bw() +
theme(legend.position=c(0.8,0.1),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Treatment") +
facet_wrap(~ genotype, ncol = 3)biomass.facet.plot = biomass.facet.plot + labs(title='Supplemental Figure S3')
print(biomass.facet.plot)
Treatment differences in biomass
Statistical analysis for fresh-weight biomass differences. Use pairwise comparisons on S. viridis and S. italica over two-day intervals.
analyze_biomass = function(vis.data, genotype) {
days = c()
diff.low = c()
diff.up = c()
pvals = c()
for(day in levels(factor(as.integer(vis.data$dap)))) {
day = as.integer(day)
control = vis.data[(as.integer(vis.data$dap) == day |
as.integer(vis.data$dap) == day + 1) &
vis.data$genotype == genotype &
vis.data$treatment == 100,]$fw_biomass
drought = vis.data[(as.integer(vis.data$dap) == day |
as.integer(vis.data$dap) == day + 1) &
vis.data$genotype == genotype &
vis.data$treatment == 33,]$fw_biomass
# If the control group has a lot more replicates (e.g. A10 and B100), randomly sample to drought sample size
if (genotype == 'A10' | genotype == 'B100') {
control = sample(control, size = length(drought))
}
test = t.test(x=control, y=drought)
days = c(days, day)
diff.low = c(diff.low, test$conf.int[1])
diff.up = c(diff.up, test$conf.int[2])
pvals = c(pvals, test$p.value)
}
results = data.frame(dap=as.numeric(days),
conf.int.low=diff.low,
conf.int.up=diff.up,
pvalue=pvals)
return(results)
}Test for treatment effect on two-day intervals
a10.biomass.results = analyze_biomass(vis.data, 'A10')
b100.biomass.results = analyze_biomass(vis.data, 'B100')Control for multiple testing by controlling the FDR
qvalues.biomass = qvalue(c(a10.biomass.results$pvalue, b100.biomass.results$pvalue))
a10.biomass.results$qvalue = qvalues.biomass$qvalues[1:nrow(a10.biomass.results)]
b100.biomass.results$qvalue = qvalues.biomass$qvalues[
(nrow(a10.biomass.results)+1):(nrow(a10.biomass.results) +
nrow(b100.biomass.results))]Assign genotypes for merged table
a10.biomass.results$genotype = 'A10'
b100.biomass.results$genotype = 'B100'
print(a10.biomass.results)## dap conf.int.low conf.int.up pvalue qvalue genotype
## 1 11 0.02914367 0.2201204 1.323444e-02 1.515363e-02 A10
## 2 12 -0.14788575 0.1756949 8.610677e-01 5.049914e-01 A10
## 3 13 -0.16165563 0.2625467 6.306207e-01 3.790869e-01 A10
## 4 14 -0.18405590 0.4646706 3.796492e-01 2.552294e-01 A10
## 5 15 -0.27354779 0.6885945 3.821226e-01 2.552294e-01 A10
## 6 16 -0.39537050 0.6661931 5.981374e-01 3.784843e-01 A10
## 7 17 0.10518290 3.1083996 3.780833e-02 3.952667e-02 A10
## 8 18 -0.13454475 2.2779082 7.835338e-02 7.138478e-02 A10
## 9 19 0.55639946 3.2418022 7.520995e-03 9.518142e-03 A10
## 10 20 1.79793954 4.9902898 2.723816e-04 4.093438e-04 A10
## 11 21 2.01188934 6.2302349 6.230113e-04 8.812062e-04 A10
## 12 22 4.48502027 8.5942426 2.003136e-06 4.013837e-06 A10
## 13 23 4.38451816 10.3131422 2.161645e-04 3.465163e-04 A10
## 14 25 6.34140762 10.9450081 1.001322e-07 3.009637e-07 A10
## 15 26 6.34064811 11.1993657 2.033517e-07 5.432952e-07 A10
## 16 27 8.43389988 12.1214602 4.849913e-11 2.332354e-10 A10
## 17 28 9.48931032 12.8236353 6.467846e-13 7.776071e-12 A10
## 18 29 8.43000027 11.5285952 1.588593e-12 1.041987e-11 A10
## 19 30 10.37389366 13.1633241 1.396552e-14 3.358054e-13 A10
## 20 31 9.37549132 12.8536239 1.733372e-12 1.041987e-11 A10
## 21 32 8.77790792 12.5759125 8.718795e-11 3.494103e-10 A10
## 22 33 9.13675701 14.3597517 1.794587e-06 4.013837e-06 A10print(b100.biomass.results)## dap conf.int.low conf.int.up pvalue qvalue genotype
## 1 11 -0.1298704 -0.02014007 9.958668e-03 1.197297e-02 B100
## 2 12 -0.0805700 0.14440699 5.635119e-01 3.662114e-01 B100
## 3 13 -0.2655002 0.05797019 1.994123e-01 1.453010e-01 B100
## 4 14 -0.2982824 0.27327839 9.285308e-01 5.315909e-01 B100
## 5 15 -0.5579089 0.03474439 8.015652e-02 7.138478e-02 B100
## 6 16 -0.4525040 0.74118118 6.175631e-01 3.790869e-01 B100
## 7 17 -1.0828481 1.08346063 9.995283e-01 5.462267e-01 B100
## 8 18 -0.2717842 1.41231582 1.737604e-01 1.347782e-01 B100
## 9 19 -1.1177500 1.15611705 9.725059e-01 5.438189e-01 B100
## 10 20 -0.2312011 2.65901500 9.508123e-02 8.165209e-02 B100
## 11 21 -0.2643127 2.63411000 1.034254e-01 8.575509e-02 B100
## 12 22 -0.5396748 3.20099463 1.535486e-01 1.230708e-01 B100
## 13 23 -2.9434030 7.26769789 3.233459e-01 2.286751e-01 B100
## 14 25 -1.2656478 5.61185394 1.982686e-01 1.453010e-01 B100
## 15 26 -0.1065401 5.98201981 5.765873e-02 5.545690e-02 B100
## 16 27 0.1611022 6.85503580 4.110506e-02 4.118268e-02 B100
## 17 28 0.7754796 9.21814913 2.356266e-02 2.575326e-02 B100
## 18 29 4.0254297 10.07486231 1.391536e-04 2.389995e-04 B100
## 19 30 6.9552625 11.09939185 7.806236e-09 2.681477e-08 B100
## 20 31 5.5617792 10.73510877 2.262272e-06 4.184388e-06 B100
## 21 32 3.0137106 11.18548970 2.080233e-03 2.778882e-03 B100
## 22 33 8.1298505 13.84599687 1.865143e-06 4.013837e-06 B100Output the A10 and B100 tables to the same file
write.table(a10.biomass.results, file='biomass.stats.csv', sep = ',',
row.names = FALSE, append = FALSE)
write.table(b100.biomass.results, file='biomass.stats.csv', sep = ',',
row.names = FALSE, append = TRUE, col.names = FALSE)Biomass response to drought
Calculate the difference in biomass per genotype per day between the 33% and 100% water treatment groups.
drought_response = function(vis.data, genotypes) {
# Initialize drought response data frame with days
drought_resp = data.frame(day=c(
min(as.integer(vis.data$dap)):max(as.integer(vis.data$dap))))
# Initialize genotypes
for(g in genotypes) {
control = paste(g,'control',sep='')
drought = paste(g,'drought',sep='')
# Calculate median biomass per treatment per day
drought_resp[,paste(control)] = 0
drought_resp[,paste(drought)] = 0
for(d in min(as.integer(vis.data$dap)):max(as.integer(vis.data$dap))) {
drought_resp[drought_resp$day == d, paste(control)] =
median(vis.data[vis.data$genotype == g &
vis.data$treatment == 100 &
as.integer(vis.data$dap) == d,]$fw_biomass)
drought_resp[drought_resp$day == d, paste(drought)] =
median(vis.data[vis.data$genotype == g &
vis.data$treatment == 33 &
as.integer(vis.data$dap) == d,]$fw_biomass)
}
}
drought_resp = drought_resp[-14,]
# Calculate biomass loss to drought
days = c()
response = c()
genotype=c()
for(r in 1:nrow(drought_resp)) {
days = c(days, rep(drought_resp[r,]$day,length(genotypes)))
for(g in genotypes) {
control = paste(g,'control',sep='')
drought = paste(g,'drought',sep='')
response = c(response, drought_resp[r,paste(drought)] -
drought_resp[r,paste(control)])
genotype = c(genotype, g)
}
}
dr.df = data.frame(day=days,response=response,genotype=as.factor(genotype))
return(dr.df)
}A10 vs B100 drought responses
drought.response.parents = drought_response(vis.data,c('A10','B100'))Plot drought response results
genotype.colors = c('#E6A024','#4AB859')
resp.plot = ggplot(drought.response.parents,
aes(x=day, y=response, group=genotype, color=genotype)) +
geom_point(size=2.5) +
geom_smooth(method="loess") +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(name="Reduced accumulated biomass (g)") +
theme_bw() +
theme(legend.position=c(0.8,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
scale_colour_manual(values=genotype.colors) +
labs(color='Genotype')resp.plot = resp.plot + labs(title='Figure 4C')
print(resp.plot)
Growth curve analysis
In this experiment plants appear to follow relatively basic asymptotic growth patterns. Here we use non-linear least squares regression analysis to model a three-component logistic growth function for each genotype x treatment group.
First load three functions for non-linear growth curve analysis reproduced from:
Paine CET, Marthews TR, Vogt DR, Purves D, Rees M, Hector A, Turnbull LA (2012) How to fit nonlinear plant growth models and calculate growth rates: an update for ecologists. Methods in Ecology and Evolution 3: 245–256. doi: 10.1111/j.2041-210X.2011.00155.x.
output.logis.nlsList <- function(fit, times, CI = F, LOG = F, alpha = 0.05){
coef <- coef(fit)
params <- transform_param.logis(coef)
rates <- list()
groups <- rownames(params)
n.groups <- nrow(coef)
# compute rates for each group seperately
rates <- list()
for(i in 1:(n.groups)){
K <- params[i,1]; r <- params[i,2]; M0 <- params[i,3]
rates[[i]] = data.frame(
times = times,
M = (M0*K)/(M0+(K-M0)*exp(-r*times)),
AGR = (r*M0*K*(K-M0)*exp(-r*times))/(M0+(K-M0)*exp(-r*times))^2
)
rates[[i]]$RGRt <- rates[[i]]$AGR/rates[[i]]$M
rates[[i]]$RGRm <- r*(1 - rates[[i]]$M/K)
if(LOG == T){
rates[[i]]$RGRt <- rates[[i]]$AGR
rates[[i]]$RGRm <- r*rates[[i]]$M*(1-rates[[i]]$M/K)
rates[[i]]$AGR <- rates[[i]]$AGR*exp(rates[[i]]$M)
}
# commute CIs for each group's estaimates, if desired
if(CI == T){
cov <- summary(fit)$cov[i,,]
x <- y <- data.frame(rmvnorm(n=1000, mean=c(coef[i, 1], coef[i, 2], coef[i, 3]), sigma=cov))
x$K <- y[,1]
x$r <- 1/y[,3]
x$M0 <- y[,1]/(1 + exp(y[,2]/y[,3]))
M <- AGR <- RGRt <- RGRm <- matrix(NA, ncol = length(times), nrow = nrow(x))
for(j in 1:nrow(x)){
K <- x[j,4]; r <- x[j,5]; M0 <- x[j,6]
M[j,] <- (M0*K)/(M0+(K-M0)*exp(-r*times))
AGR[j,] <- (r*M0*K*(K-M0)*exp(-r*times))/(M0+(K-M0)*exp(-r*times))^2
RGRt[j,] <- AGR[j,]/M[j,]
RGRm[j,] <- r*(1 - M[j,]/K)
if(LOG ==T){
RGRt[j,] <- AGR[j,]
RGRm[j,] <- r*M[j,]*(1 - M[j,]/K)
AGR[j,] <- AGR[j,]*exp(M[j,])
}
}
}
CIs <- summarizer(list(M = M, AGR = AGR, RGRt = RGRt, RGRm = RGRm), alpha)
rates[[i]] <- cbind(rates[[i]], CIs)
}
names(rates) <- rownames(params)
# now compute differences among groups
diffs <- list()
for(i in 1:(n.groups-1)){
Ki <- params[i,1]; ri <- params[i,2]; M0i <- params[i,3]
for(j in (i+1):n.groups){
Kj <- params[j,1]; rj <- params[j,2]; M0j <- params[j,3]
diffs.ij = data.frame(
times = times,
diffM = rates[[i]]$M - rates[[j]]$M,
diffAGR = rates[[i]]$AGR - rates[[j]]$AGR,
diffRGRt = rates[[i]]$RGRt - rates[[j]]$RGRt
)
# comparing RGRm has to be done on a biomass basis. So it needs special treatment. First, we aev to know what range of biomasses are shared between groups.
Mmin <- max(min(rates[[i]]$M), min(rates[[j]]$M)) # yieds the range of overlapping masses between groups i & j
Mmax <- min(max(rates[[i]]$M), max(rates[[j]]$M))
diffs.ij$Mseq <- seq(Mmin, Mmax, length = 100)
if(LOG == F){
diffs.ij$diffRGRm <- ri*(1 - diffs.ij$Mseq/Ki) - rj*(1 - diffs.ij$Mseq/Kj)
} else{
diffs.ij$diffRGRm <- ri*diffs.ij$Mseq*(1 - diffs.ij$Mseq/Ki) - rj*Mseq*(1 - diffs.ij$Mseq/Kj)
}
}
if(CI == T){
# get params for group i
covi <- summary(fit)$cov[i,,]
xi <- yi <- data.frame(rmvnorm(n=1000, mean=c(coef[i, 1], coef[i, 2], coef[i, 3]), sigma=covi))
xi$K <- yi[,1]
xi$r <- 1/yi[,3]
xi$M0 <- yi[,1]/(1 + exp(yi[,2]/yi[,3]))
# get params for group j
covj <- summary(fit)$cov[j,,]
xj <- yj <- data.frame(rmvnorm(n=1000, mean=c(coef[j, 1], coef[j, 2], coef[j, 3]), sigma=covj))
xj$K <- yj[,1]
xj$r <- 1/yj[,3]
xj$M0 <- yj[,1]/(1 + exp(yj[,2]/yj[,3]))
# now compute diffs for each random set of drawn parameters
Mi <- Mj <- AGRi <- AGRj <- RGRti <- RGRtj <- RGRmi <- RGRmj <- diffM <- diffAGR <- diffRGRt <- diffRGRm <- matrix(NA, ncol = length(times), nrow = nrow(xi))
for(k in 1:nrow(xi)){
Ki <- xi[k,4]; ri <- xi[k,5]; M0i <- xi[k,6]
Kj <- xj[k,4]; rj <- xj[k,5]; M0j <- xj[k,6]
Mi[k,] <- (M0i*Ki)/(M0i+(Ki-M0i)*exp(-ri*times))
Mj[k,] <- (M0j*Kj)/(M0j+(Kj-M0j)*exp(-rj*times))
AGRi[k,] <- (ri*M0i*Ki*(Ki-M0i)*exp(-ri*times))/(M0i+(Ki-M0i)*exp(-ri*times))^2
AGRj[k,] <- (rj*M0j*Kj*(Kj-M0j)*exp(-rj*times))/(M0j+(Kj-M0j)*exp(-rj*times))^2
RGRti[k,] <- AGRi[k,]/Mi[k,]
RGRtj[k,] <- AGRj[k,]/Mj[k,]
RGRmi[k,] <- ri*(1 - diffs.ij$Mseq/Ki)
RGRmj[k,] <- rj*(1 - diffs.ij$Mseq/Kj)
if(LOG == T){
RGRti[k,] <- AGRi[k,]
RGRtj[k,] <- AGRj[k,]
RGRmi[k,] <- ri*diffs.ij$Mseq*(1 - diffs.ij$Mseq/Ki)
RGRmj[k,] <- rj*diffs.ij$Mseq*(1 - diffs.ij$Mseq/Kj)
AGRi[k,] <- AGRi[k,]*exp(Mi[k,])
AGRj[k,] <- AGRj[k,]*exp(Mj[k,])
}
diffM[k,] <- Mi[k,] - Mj[k,]
diffAGR[k,] <- AGRi[k,] - AGRj[k,]
diffRGRt[k,] <- RGRti[k,] - RGRtj[k,]
diffRGRm[k,] <- RGRmi[k,] - RGRmj[k,]
}
CIs <- summarizer(list(diffM = diffM, diffAGR = diffAGR, diffRGRt = diffRGRt, diffRGRm = diffRGRm), alpha)
diffs[[paste(groups[i], groups[j], sep = "_")]] <- cbind(diffs.ij, CIs)
} else{
diffs[[paste(groups[i], groups[j], sep = "_")]] <- diffs.ij
}
} # end loop over pairwise combinations of groups
out <- list(params = params, rates = rates, diffs = diffs)
return(out)
}
transform_param.logis <- function(coef){
K = coef[1]
r = 1/(coef[3])
M0 = K/(1 + exp(coef[2]/coef[3])) #untransform best-fit parameters to K, r and M0
if(is.data.frame(K)){
out <- cbind(K, r, M0)
} else {
out <- c(K, r, M0)
}
names(out) <- c("K", "r", "M0")
return(out)
}
# this function returns confidence envelopes around growth trajectories, and growth rates.
summarizer <- function(dat, alpha){
n <- length(dat)
quantiles <- c(alpha/2, 1-(alpha/2))
CIs <- data.frame(matrix(NA, ncol(dat[[1]]), n*2))
names(CIs) <- paste(rep(names(dat), each = 2), c("lo", "hi"), sep = ".")
for(i in 1:n){
CIs[,(2*i-1):(2*i)] <- t(apply(dat[[i]], 2, quantile, quantiles, na.rm = T))
}
return(CIs)
}Subset the VIS data
sub.vis.data = vis.data[(vis.data$treatment == 100 | vis.data$treatment == 33),
c("plant_id","dap","fw_biomass","group")]
sub.vis.data$group = factor(sub.vis.data$group)Group data
grouped.sub.vis.data = groupedData(fw_biomass ~ dap | group, sub.vis.data)Fit three-component logistic functions for each genotype x treatment group
fit.logis = nlsList(fw_biomass ~ SSlogis(dap, Asym, xmid, scal),
data = grouped.sub.vis.data)Output estimated growth rate parameters at even time intervals
est.interval = seq(min(sub.vis.data$dap), max(sub.vis.data$dap), length = 100)
out.fit.logis = output.logis.nlsList(fit.logis, times = est.interval,
CI = TRUE, LOG = FALSE, alpha = 0.05)Compute the time and magnitude of maximum growth rate
logis.results = data.frame(group=levels(sub.vis.data$group))
logis.results$max.AGR = 0
logis.results$max.AGR.dap = 0
group.order = names(out.fit.logis$rates)
for(i in 1:length(group.order)) {
logis.results[logis.results == group.order[i],]$max.AGR.dap =
out.fit.logis$rates[[i]]$times[out.fit.logis$rates[[i]]$AGR ==
max(out.fit.logis$rates[[i]]$AGR)]
logis.results[logis.results == group.order[i],]$max.AGR =
max(out.fit.logis$rates[[i]]$AGR)
}
print(logis.results)## group max.AGR max.AGR.dap
## 1 A10-100 2.994229 24.23556
## 2 A10-33 1.737080 25.36054
## 3 B100-100 2.389559 26.26053
## 4 B100-33 1.631987 24.68555
## 5 R102-100 2.593812 28.06050
## 6 R102-33 1.826986 26.48553
## 7 R128-100 2.995409 25.36054
## 8 R128-33 1.781943 24.23556
## 9 R133-100 3.267807 23.11057
## 10 R133-33 2.207510 23.33557
## 11 R161-100 2.858742 25.58554
## 12 R161-33 1.750648 24.68555
## 13 R187-100 2.444363 24.46055
## 14 R187-33 1.617136 24.46055
## 15 R20-100 3.240594 22.88557
## 16 R20-33 1.869087 22.66058
## 17 R70-100 2.878034 23.78556
## 18 R70-33 1.592505 23.33557
## 19 R98-100 2.509740 24.23556
## 20 R98-33 1.481989 24.68555Plot logistic growth models for S. viridis and S. italica and absolute growth rates over time
parent.groups = c("A10-100", "A10-33", "B100-100", "B100-33")
group.colors = c("#F8766D", "#7CAE00", "#00BFC4", "#C77CFF")
# Biomass over time
gca.plot = ggplot(vis.data[(vis.data$genotype == 'A10' | vis.data$genotype == 'B100') &
(vis.data$treatment == 100 | vis.data$treatment == 33),],
aes(x=dap, y=fw_biomass, color=factor(group))) +
geom_point(size=2.5) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(-1,41), name="Estimated biomass (g)") +
theme_bw() +
theme(legend.position=c(0.2,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
scale_colour_manual(name='Genotype-treatment',
values=group.colors, labels=parent.groups)
agr.plot = ggplot() +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(name="Absolute growth rate (g/day)") +
theme_bw() +
theme(legend.position=c(0.2,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
scale_colour_manual(name='Genotype-treatment',
values=group.colors, labels=parent.groups)
# Add logistic growth models and confidence intervals
for(g in 1:length(parent.groups)) {
group = parent.groups[g]
i = match(group, group.order)
group.rates = as.data.frame(out.fit.logis$rates[i])
col.name = gsub('-', '.', group)
conf.int = data.frame(x=c(group.rates[,paste(col.name, '.times', sep='')],
rev(group.rates[,paste(col.name, '.times', sep='')])),
y=c(group.rates[,paste(col.name, '.M.lo', sep='')],
rev(group.rates[,paste(col.name, '.M.hi', sep='')])))
agr.conf.int = data.frame(x=c(group.rates[,paste(col.name, '.times', sep='')],
rev(group.rates[,paste(col.name, '.times', sep='')])),
y=c(group.rates[,paste(col.name, '.AGR.lo', sep='')],
rev(group.rates[,paste(col.name, '.AGR.hi', sep='')])))
gca.plot = gca.plot + geom_polygon(data=conf.int, aes(x=x, y=y),
fill='gray60', color=NA, alpha=0.4) +
geom_line(data=group.rates, aes_string(x=paste(col.name, '.times', sep=''),
y=paste(col.name, '.M', sep='')),
color=group.colors[g])
agr.plot = agr.plot +
geom_polygon(data=agr.conf.int, aes(x=x, y=y),
fill='gray60', color=NA, alpha=0.4) +
geom_line(data=group.rates, aes_string(x=paste(col.name, '.times', sep=''),
y=paste(col.name, '.AGR', sep='')),
color=group.colors[g])
}gca.plot = gca.plot + labs(title='Figure 4B')
print(gca.plot)
agr.plot = agr.plot + labs(title='Figure 4C')
print(agr.plot)
Plot logistic growth models for all Setaria genotypes
plot_gca = function(genotype) {
groups = c(paste(genotype,'-100',sep=''), paste(genotype,'-33',sep=''))
group.colors = c("#00BFC4", "#F8766D")
# Biomass over time
gca.plot = ggplot(vis.data[(vis.data$genotype == genotype) &
(vis.data$treatment == 100 | vis.data$treatment == 33),],
aes(x=dap, y=fw_biomass, color=factor(treatment))) +
geom_point(size=1) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(-1,41), name="Estimated biomass (g)") +
theme_bw() +
theme(legend.position='none',
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Treatment") +
facet_grid(. ~ genotype)
# Add logistic growth models and confidence intervals
for(g in 1:length(groups)) {
group = groups[g]
i = match(group, group.order)
group.rates = as.data.frame(out.fit.logis$rates[i])
col.name = gsub('-', '.', group)
conf.int = data.frame(x=c(group.rates[,paste(col.name, '.times', sep='')],
rev(group.rates[,paste(col.name, '.times', sep='')])),
y=c(group.rates[,paste(col.name, '.M.lo', sep='')],
rev(group.rates[,paste(col.name, '.M.hi', sep='')])))
gca.plot = gca.plot + geom_polygon(data=conf.int, aes(x=x, y=y),
fill='gray60', color=NA, alpha=0.4) +
geom_line(data=group.rates, aes_string(x=paste(col.name, '.times', sep=''),
y=paste(col.name, '.M', sep='')),
color=group.colors[g])
}
return(gca.plot)
}Multiple plot function modified from the Cookbook for R by Winston Chang: http://www.cookbook-r.com/Graphs/Multiple_graphs_on_one_page_(ggplot2)/
multiplot <- function(..., plotlist=NULL, cols=1, layout=NULL) {
# Make a list from the ... arguments and plotlist
plots <- c(list(...), plotlist)
numPlots = length(plots)
# If layout is NULL, then use 'cols' to determine layout
if (is.null(layout)) {
# Make the panel
# ncol: Number of columns of plots
# nrow: Number of rows needed, calculated from # of cols
layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
ncol = cols, nrow = ceiling(numPlots/cols),byrow=TRUE)
}
if (numPlots==1) {
print(plots[[1]])
} else {
# Set up the page
grid.newpage()
pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
# Make each plot, in the correct location
for (i in 1:numPlots) {
# Get the i,j matrix positions of the regions that contain this subplot
matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
layout.pos.col = matchidx$col))
}
}
}Generate the Setaria biomass/growth plots
a10.gca.plot = plot_gca('A10')
b100.gca.plot = plot_gca('B100')
r102.gca.plot = plot_gca('R102')
r128.gca.plot = plot_gca('R128')
r133.gca.plot = plot_gca('R133')
r161.gca.plot = plot_gca('R161')
r187.gca.plot = plot_gca('R187')
r20.gca.plot = plot_gca('R20')
r70.gca.plot = plot_gca('R70')
r98.gca.plot = plot_gca('R98')Supplemental Figure S3
multiplot(a10.gca.plot, b100.gca.plot, r102.gca.plot, r128.gca.plot,
r133.gca.plot, r161.gca.plot, r187.gca.plot, r20.gca.plot,
r70.gca.plot, r98.gca.plot, cols=3)
agr.plot = agr.plot + labs(title='Figure 4C')
print(agr.plot)
Water use efficiency
In this section we combine data on the volume of water added to each plant per water application (up to the target weight) with fresh-weight biomass data to estimate water use efficiency.
Download water volume data (n = 33,496 water applications).
if (!file.exists('B2_watering_data_phenofront.csv')) {
download.file('http://files.figshare.com/1845363/B2_watering_data_phenofront.csv',
'B2_watering_data_phenofront.csv')
}Read data and format for analysis
water.data = read.table(file="B2_watering_data_phenofront.csv", sep=",", header=TRUE)Remove the empty pot controls (Barcodes start with zero).
water.data = water.data[grep('000A', water.data$plant.barcode, invert = TRUE),]Remove snapshots that were not valid waterings (water.amount = -1).
water.data = water.data[water.data$water.amount != -1,]Add a treatment group column coded in the barcode.
water.data$treatment <- NA
water.data$treatment[grep("AA", water.data$plant.barcode)] <- 100
water.data$treatment[grep("AB", water.data$plant.barcode)] <- 0
water.data$treatment[grep("AC", water.data$plant.barcode)] <- 16
water.data$treatment[grep("AD", water.data$plant.barcode)] <- 33
water.data$treatment[grep("AE", water.data$plant.barcode)] <- 66Add a plant genotype column coded in the barcode.
water.data$genotype <- NA
water.data$genotype[grep("p1", water.data$plant.barcode)] <- 'A10'
water.data$genotype[grep("p2", water.data$plant.barcode)] <- 'B100'
water.data$genotype[grep("r1", water.data$plant.barcode)] <- 'R20'
water.data$genotype[grep("r2", water.data$plant.barcode)] <- 'R70'
water.data$genotype[grep("r3", water.data$plant.barcode)] <- 'R98'
water.data$genotype[grep("r4", water.data$plant.barcode)] <- 'R102'
water.data$genotype[grep("r5", water.data$plant.barcode)] <- 'R128'
water.data$genotype[grep("r6", water.data$plant.barcode)] <- 'R133'
water.data$genotype[grep("r7", water.data$plant.barcode)] <- 'R161'
water.data$genotype[grep("r8", water.data$plant.barcode)] <- 'R187'Add a genotype x treatment group column.
water.data$group = paste(water.data$genotype,'-',water.data$treatment,sep='')Encode the calendar time column as a date-time.
water.data$timestamp = ymd_hms(water.data$timestamp)Add a column for days after planting since the planting date.
water.data$dap = as.numeric(water.data$timestamp - planting_date)Some plants were removed from the system but were not inactivated, so remove them.
bad.cars = unique(water.data$car.tag[water.data$water.amount >120 &
water.data$dap > 14])
water.data = water.data[!water.data$car.tag %in% bad.cars,]Data for Figure 1B
Extract water data for S. viridis for Figure 1B.
sv.water = water.data[(water.data$treatment == 100 | water.data$treatment == 33) &
water.data$genotype == 'A10',]Classify water job as early or later in the day. Early water treatments started at ZT23, later treatments started at ZT8.
sv.water$early.late <- NA #restrict to the scheduled watering later in the run
sv.water$early.late[hour(sv.water$timestamp) < 7] <- 'ZT23'
sv.water$early.late[hour(sv.water$timestamp) < 15 &
hour(sv.water$timestamp) > 11 ] <- 'ZT8'
sv.water = sv.water[!is.na(sv.water$early.late),]
sv.water$dap = day(sv.water$timestamp) + 4
sv.water = sv.water[sv.water$dap > 14,]Plot the S. viridis water volume data
water.plot = ggplot(sv.water, aes(y = water.amount, x = dap,
group = factor(interaction(early.late,
treatment,dap)),
fill = factor(interaction(early.late,
treatment)),
color = factor(interaction(early.late,
treatment)))) +
geom_boxplot(outlier.colour = NULL) +
scale_x_continuous("Days after planting") +
scale_y_continuous("Water volume (ml)") +
theme_bw() +
theme(legend.position = c(0.2,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(fill='Water treatment', color='Water treatment')water.plot = water.plot + labs(title='Figure 1B')
print(water.plot)
WUE function
Calculate WUE for S. viridis and S. italica
wue.parents = function(water, biomass) {
# WUE data vectors
dap.list = c()
plant.list = c()
water.list = c()
biomass.list = c()
treatment.list = c()
group.list=c()
# Get unique barcodes for biomass
barcodes = unique(biomass[(biomass$genotype=='A10' | biomass$genotype=='B100') &
(biomass$treatment == 100 |
biomass$treatment == 33),]$plant_id)
for(barcode in barcodes) {
snapshots = biomass[biomass$plant_id==barcode,]
snapshots = snapshots[with(snapshots, order(dap)),]
for(row in 1:nrow(snapshots)) {
total.water = sum(water[water$dap <= snapshots[row,]$dap &
water$plant.barcode == barcode,]$water.amount)
dap.list = c(dap.list, snapshots[row,]$dap)
plant.list = c(plant.list, barcode)
water.list = c(water.list, total.water)
biomass.list = c(biomass.list, snapshots[row,]$fw_biomass)
treatment.list = c(treatment.list, snapshots[row,]$treatment)
group.list = c(group.list,snapshots[row,]$group)
}
}
wue.data = data.frame(plant_id=plant.list,
dap=dap.list,
water=water.list,
biomass=biomass.list,
treatment=treatment.list,
group=group.list)
return(wue.data)
}
parents.wue = wue.parents(water.data, vis.data)
parents.wue = parents.wue[parents.wue$water != 0,]Plot S. viridis and S. italica WUE in mg/mL.
wue.plot = ggplot(parents.wue, aes(x=dap, y=(biomass/water)*1000,
color=factor(group))) +
geom_point(size=2.5) +
geom_smooth(method="loess", size=1) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(-2.1, 26),
name="Water-use efficiency (mg/mL)") +
theme_bw() +
theme(legend.position=c(0.2,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color='Genotype-treatment')wue.plot = wue.plot + labs(title='Figure 4D')
print(wue.plot)
Statistical analysis of WUE
Genotype WUE function
wue = function(water, biomass, genotype) {
# WUE data vectors
dap.list = c()
plant.list = c()
water.list = c()
biomass.list = c()
treatment.list = c()
# Get unique barcodes for biomass
barcodes = unique(biomass[biomass$genotype==genotype &
(biomass$treatment == 100 |
biomass$treatment == 33),]$plant_id)
for(barcode in barcodes) {
snapshots = biomass[biomass$plant_id==barcode,]
snapshots = snapshots[with(snapshots, order(dap)),]
for(row in 1:nrow(snapshots)) {
total.water = sum(water[water$dap <= snapshots[row,]$dap &
water$plant.barcode == barcode,]$water.amount)
dap.list = c(dap.list, snapshots[row,]$dap)
plant.list = c(plant.list, barcode)
water.list = c(water.list, total.water)
biomass.list = c(biomass.list, snapshots[row,]$fw_biomass)
treatment.list = c(treatment.list, snapshots[row,]$treatment)
}
}
wue.data = data.frame(plant_id=plant.list,
dap=dap.list,
water=water.list,
biomass=biomass.list,
treatment=treatment.list)
return(wue.data)
}Analyse pair-wise treatment differences for two-day increments for each genotype.
analyze_wue = function(wue.data) {
days = c()
diff.low = c()
diff.up = c()
pvals = c()
for(day in levels(factor(as.integer(wue.data$dap)))) {
day = as.integer(day)
control = wue.data[(as.integer(wue.data$dap) == day |
as.integer(wue.data$dap) == day + 1) &
wue.data$treatment == 100,]
drought = wue.data[(as.integer(wue.data$dap) == day |
as.integer(wue.data$dap) == day + 1) &
wue.data$treatment == 33,]
control.wue = control$biomass / control$water
drought.wue = drought$biomass / drought$water
test = t.test(x=control.wue,y=drought.wue)
days = c(days,day)
diff.low = c(diff.low, test$conf.int[1])
diff.up = c(diff.up, test$conf.int[2])
pvals = c(pvals, test$p.value)
}
results = data.frame(dap=as.numeric(days),
conf.int.low=diff.low,
conf.int.up=diff.up,
pvalue=pvals)
return(results)
}Analyze WUE for each genotype and aggregate results
wue.results = data.frame()
for(genotype in c('A10', 'B100')) {
genotype.wue = wue(water.data, vis.data, genotype)
genotype.wue = genotype.wue[genotype.wue$water != 0,]
genotype.wue.results = analyze_wue(genotype.wue)
genotype.wue.results$genotype = genotype
wue.results = rbind(wue.results, genotype.wue.results)
}Control for multiple testing by controlling the FDR
qvalues.wue = qvalue(wue.results$pvalue)
wue.results$qvalue = qvalues.wue$qvalues
write.table(wue.results, file='wue.stats.csv', sep=',', row.names=FALSE)
print(wue.results)## dap conf.int.low conf.int.up pvalue genotype qvalue
## 1 11 0.0001290017 6.743361e-04 0.0055453872 A10 0.04066617
## 2 12 -0.0003304552 6.531982e-04 0.5007434856 A10 0.90056523
## 3 13 -0.0004887093 6.645100e-04 0.7558400553 A10 0.95533507
## 4 14 -0.0005521785 8.912840e-04 0.6285589728 A10 0.92188649
## 5 15 -0.0006900247 1.272981e-03 0.5413573872 A10 0.91614327
## 6 16 -0.0013750018 1.352784e-03 0.9863361510 A10 0.98633615
## 7 17 -0.0032500144 1.960050e-03 0.5946773334 A10 0.92188649
## 8 18 -0.0026416988 1.024200e-03 0.3665557006 A10 0.80642254
## 9 19 -0.0025138146 1.198806e-03 0.4687644330 A10 0.90056523
## 10 20 -0.0021111216 2.216305e-03 0.9594528855 A10 0.98633615
## 11 21 -0.0013777865 3.004615e-03 0.4446854076 A10 0.90056523
## 12 22 -0.0012554051 3.263703e-03 0.3590888551 A10 0.80642254
## 13 23 -0.0029737246 4.659146e-03 0.6177633347 A10 0.92188649
## 14 25 -0.0016604803 2.491094e-03 0.6763471555 A10 0.95533507
## 15 26 -0.0016183794 2.333780e-03 0.7068814380 A10 0.95533507
## 16 27 -0.0014579106 1.859317e-03 0.8033499452 A10 0.95533507
## 17 28 -0.0015730725 1.687446e-03 0.9422861297 A10 0.98633615
## 18 29 -0.0016425521 1.397919e-03 0.8687331415 A10 0.98633615
## 19 30 -0.0017235963 1.223038e-03 0.7270761809 A10 0.95533507
## 20 31 -0.0017633197 1.027390e-03 0.5898632801 A10 0.92188649
## 21 32 -0.0018009970 9.249981e-04 0.5116847871 A10 0.90056523
## 22 33 -0.0019234072 2.482920e-03 0.7841395855 A10 0.95533507
## 23 11 -0.0004663250 2.293599e-04 0.4748570604 B100 0.90056523
## 24 12 -0.0003827145 4.190170e-04 0.9249400490 B100 0.98633615
## 25 13 -0.0007208859 1.779908e-04 0.2235222936 B100 0.64019648
## 26 14 -0.0007089761 6.505828e-04 0.9282507680 B100 0.98633615
## 27 15 -0.0013128372 2.973716e-04 0.1995779600 B100 0.62724502
## 28 16 -0.0014265056 1.455777e-03 0.9826344845 B100 0.98633615
## 29 17 -0.0054259745 -3.416845e-04 0.0302482577 B100 0.12099303
## 30 18 -0.0039584923 7.814527e-05 0.0582767730 B100 0.21368150
## 31 19 -0.0053958070 -8.822144e-04 0.0097441669 B100 0.06124905
## 32 20 -0.0049249253 -5.837513e-04 0.0161879273 B100 0.07914098
## 33 21 -0.0062171326 -1.283934e-03 0.0055344094 B100 0.04066617
## 34 22 -0.0056491252 -1.224773e-03 0.0043669943 B100 0.04066617
## 35 23 -0.0061377032 -1.581340e-03 0.0036480639 B100 0.04066617
## 36 25 -0.0056519496 -1.646284e-03 0.0009985777 B100 0.03550789
## 37 26 -0.0050694636 -1.330337e-03 0.0016139948 B100 0.03550789
## 38 27 -0.0041897428 -5.455619e-04 0.0130667633 B100 0.07186720
## 39 28 -0.0036907493 -2.531025e-04 0.0261291938 B100 0.11496845
## 40 29 -0.0033005146 1.619285e-04 0.0734489761 B100 0.24859653
## 41 30 -0.0026887815 6.877754e-04 0.2327987196 B100 0.64019648
## 42 31 -0.0025257356 7.093671e-04 0.2546041527 B100 0.65897545
## 43 32 -0.0024342851 8.292998e-04 0.3164971519 B100 0.77365970
## 44 33 -0.0028628094 2.203841e-03 0.7831321502 B100 0.95533507Model tiller number
In this section we use linear modeling to estimate tiller number.
Calculate height-width ratio.
vis.data$height_width_ratio = vis.data$height_above_bound / vis.data$extent_xDownload data for manually measured tiller counts 195 random images.
if (!file.exists('200tiller_counts.txt')) {
download.file('http://files.figshare.com/1845359/200tiller_counts.txt',
'200tiller_counts.txt')
}Read manual tiller count data.
tiller.counts = read.table(file="200tiller_counts.txt",sep="\t",header=TRUE)Format date.
tiller.counts$date = as.POSIXct(paste(paste(tiller.counts$year,
tiller.counts$month,
tiller.counts$day, sep='-'),
paste(tiller.counts$hours,
tiller.counts$minutes,
tiller.counts$seconds, sep=':'), sep=' '))Merge with VIS snapshot table.
tiller.counts = merge(x=tiller.counts, y=vis.data, by = 'date')Model tiller counts.
tiller.model.full = lm(ave_tillers ~ fw_biomass + height_above_bound + solidity +
extent_x + height_width_ratio, data=tiller.counts)Variance inflation factor.
vif(tiller.model.full)## fw_biomass height_above_bound solidity
## 17.338117 22.493335 1.651753
## extent_x height_width_ratio
## 14.767073 3.280364Drop height and recalculate VIF
tiller.model1 = lm(ave_tillers ~ fw_biomass + solidity + extent_x +
height_width_ratio, data=tiller.counts)
vif(tiller.model1)## fw_biomass solidity extent_x
## 10.145172 1.393258 11.362312
## height_width_ratio
## 1.855525Drop extent x.
tiller.model2 = lm(ave_tillers ~ fw_biomass + solidity + height_width_ratio,
data=tiller.counts)
vif(tiller.model2)## fw_biomass solidity height_width_ratio
## 1.027223 1.044973 1.039641summary(tiller.model2)##
## Call:
## lm(formula = ave_tillers ~ fw_biomass + solidity + height_width_ratio,
## data = tiller.counts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.7006 -1.2510 -0.0308 0.8410 5.6311
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.17731 0.54049 9.579 < 2e-16 ***
## fw_biomass 0.22888 0.01275 17.957 < 2e-16 ***
## solidity -0.07202 1.87918 -0.038 0.969
## height_width_ratio -2.16290 0.32335 -6.689 2.42e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.699 on 191 degrees of freedom
## Multiple R-squared: 0.6501, Adjusted R-squared: 0.6446
## F-statistic: 118.3 on 3 and 191 DF, p-value: < 2.2e-16Drop solidity
tiller.model3 = lm(ave_tillers ~ fw_biomass + height_width_ratio, data=tiller.counts)
summary(tiller.model3)##
## Call:
## lm(formula = ave_tillers ~ fw_biomass + height_width_ratio, data = tiller.counts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.7003 -1.2525 -0.0264 0.8392 5.6326
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.1640 0.4136 12.487 < 2e-16 ***
## fw_biomass 0.2289 0.0126 18.173 < 2e-16 ***
## height_width_ratio -2.1650 0.3177 -6.815 1.18e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.694 on 192 degrees of freedom
## Multiple R-squared: 0.6501, Adjusted R-squared: 0.6465
## F-statistic: 178.4 on 2 and 192 DF, p-value: < 2.2e-16Plot partial regressions.
avPlots(model = tiller.model3, pch=16, layout=c(2,1), main="Figure 5B")
Download data for manually measured tiller counts for 646 random images.
if (!file.exists('660tiller_counts.txt')) {
download.file('http://files.figshare.com/1845358/660tiller_counts.txt',
'660tiller_counts.txt')
}Predict tiller count for a second set of ~660 randomly selected plants.
tiller660 = read.table(file="660tiller_counts.txt", sep="\t", header=TRUE)
tiller660 = merge(x=tiller660, y=vis.data, by='datetime')
tiller660$predicted_tiller_count = predict.lm(object = tiller.model3,
newdata=tiller660)Plot manual versus predicted tiller counts.
tiller.plot = ggplot(data=tiller660, aes(x=tiller_count, y=predicted_tiller_count,
color=factor(genotype))) +
geom_point(size=2.5) +
geom_abline(intercept=0, slope=1) +
scale_x_continuous(lim=c(-1,14), "Manual tiller count") +
scale_y_continuous(lim=c(-1,14), "Predicted tiller count") +
theme_bw() +
theme(legend.position=c(0.45,0.95),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold"),
legend.direction="horizontal") +
labs(color="Genotype")tiller.plot = tiller.plot + labs(title='Figure 5C')
print(tiller.plot)
Predict tiller counts for the whole data set.
vis.data$tiller_count = predict.lm(object = tiller.model3, newdata = vis.data)Plot tiller counts for S. viridis and S. italica water treatments 100% and 33% full-capacity.
tc.par.plot = ggplot(vis.data[(vis.data$genotype == 'A10' |
vis.data$genotype == 'B100') &
(vis.data$treatment == 100 |
vis.data$treatment == 33),],
aes(x=dap, y=tiller_count, color=factor(group))) +
geom_point(size=2.5) +
geom_smooth(method="loess",size=1) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(0,12),
name="Estimated tiller count") +
theme_bw() +
theme(legend.position=c(0.25,0.75),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Genotype-treatment")## Warning: Removed 1 rows containing missing values (stat_smooth).## Warning: Removed 1 rows containing missing values (stat_smooth).## Warning: Removed 2 rows containing missing values (geom_point).tc.par.plot = tc.par.plot + labs(title='Figure 5D')
print(tc.par.plot)## Warning: Removed 1 rows containing missing values (stat_smooth).## Warning: Removed 1 rows containing missing values (stat_smooth).## Warning: Removed 2 rows containing missing values (geom_point).
Plot estimated tiller counts for all genotypes.
tiller.facet.plot = ggplot(vis.data[vis.data$treatment == 100 |
vis.data$treatment == 33,],
aes(x=dap, y=tiller_count, color=factor(treatment))) +
geom_point(size=1) +
geom_smooth(method="loess",size=1) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(0,12),
name="Estimated tiller count") +
theme_bw() +
theme(legend.position=c(0.8,0.1),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Treatment") +
facet_wrap(~ genotype, ncol = 3)## Warning: Removed 1 rows containing missing values (stat_smooth).## Warning: Removed 1 rows containing missing values (stat_smooth).## Warning: Removed 1 rows containing missing values (stat_smooth).## Warning: Removed 1 rows containing missing values (stat_smooth).## Warning: Removed 1 rows containing missing values (stat_smooth).## Warning: Removed 3 rows containing missing values (stat_smooth).## Warning: Removed 1 rows containing missing values (stat_smooth).## Warning: Removed 1 rows containing missing values (stat_smooth).## Warning: Removed 1 rows containing missing values (stat_smooth).## Warning: Removed 2 rows containing missing values (geom_point).## Warning: Removed 2 rows containing missing values (geom_point).## Warning: Removed 1 rows containing missing values (geom_point).## Warning: Removed 4 rows containing missing values (geom_point).## Warning: Removed 1 rows containing missing values (geom_point).## Warning: Removed 1 rows containing missing values (geom_point).tiller.facet.plot = tiller.facet.plot + labs(title='Supplemental Figure S4')
print(tiller.facet.plot)## Warning: Removed 1 rows containing missing values (stat_smooth).## Warning: Removed 1 rows containing missing values (stat_smooth).## Warning: Removed 1 rows containing missing values (stat_smooth).## Warning: Removed 1 rows containing missing values (stat_smooth).## Warning: Removed 1 rows containing missing values (stat_smooth).## Warning: Removed 3 rows containing missing values (stat_smooth).## Warning: Removed 1 rows containing missing values (stat_smooth).## Warning: Removed 1 rows containing missing values (stat_smooth).## Warning: Removed 1 rows containing missing values (stat_smooth).## Warning: Removed 2 rows containing missing values (geom_point).## Warning: Removed 2 rows containing missing values (geom_point).## Warning: Removed 1 rows containing missing values (geom_point).## Warning: Removed 4 rows containing missing values (geom_point).## Warning: Removed 1 rows containing missing values (geom_point).## Warning: Removed 1 rows containing missing values (geom_point).
Trait table
Output a table of all of the traits input and added here. Export all of the VIS traits (except color).
h.table = vis.data
h.table$plant_id = as.character(h.table$plant_id)
h.table$wue = NA
for(r in 1:nrow(h.table)) {
row = h.table[r,]
total.water = sum(water.data[water.data$dap <= row$dap & water.data$plant.barcode == row$plant_id,]$water.amount)
h.table[r,]$wue = row$sv_area / total.water
}
write.table(h.table,file="vis.traits.csv",quote=FALSE,sep=",",row.names=FALSE)